/* Pivot Algorithm for Polymers
 * Author: Florian Seidler, Jonathan Zopes
 * Date: 12.12.2012
 * Issue: Computational Physics WT 12/13 , Bonn University
*/
#define TWO_DIMENSIONS //Specify Dimension
#include <iostream>
#include <math.h>
using namespace std;

/*-----------------------------------------*/
/*Tree Structure and Operations on the Tree*/
/*-----------------------------------------*/

/*Tree nodes*/
struct saw_node
{
int n; //Number of Sites
saw_node *left; //Left node pointer
saw_node *right; //Right node pointer
saw_node *parent; //Parent node pointer
int q; //Symmetry operation
#ifdef TWO_DIMENSIONS
int X[2]; //End-Vector
int B[4]; //Bounding Box: infx,infy,supx,supy
#endif
#ifdef THREE_DIMENSIONS
int X[3]; //End-Vector
int B[6]; //Bounding Box: infx,infy,infz,supx,supy,supz
#endif
#ifdef FOUR_DIMENSIONS
int X[4]; //End-Vector
int B[8]; //Bounding Box: infx,infy,infz,infw,supx,supy,supz,supw
#endif
};

/*Add tree node to correct position.*/
void add_node(int j, saw_node* root)
{
    int depth,bound,k,i;
    double l=0;
    saw_node* curr = root;
    //calculate the depth of the current tree
    depth=1; k=2;
    bound=2;
    while(bound<j) {
        bound+=pow(2,k);
        depth++;
        k++;
    }
    //Search for position to add
    for(i=1;i<depth;i++) {
        if((double)(j-(bound-pow(2,depth)))/(double)pow(2,depth)<=(1.0/(double)pow(2,i)+l)) {
            curr=curr->left;
        }
        else {
            curr=curr->right;
            l+=1/(double)pow(2,i);
        }
    }
    //Generate new node and set trivial initial parameters
    saw_node *toadd=new saw_node;
    toadd->left=NULL;
    toadd->right=NULL;
    toadd->parent=curr;
    //Initialize parameters here
    toadd->n=0;
    toadd->q=0;
    toadd->X[0]=0;
    toadd->X[1]=0;
    #ifdef THREE_DIMENSIONS
    toadd->X[2]=0;
    #endif
    #ifdef FOUR_DIMENSIONS
    toadd->X[2]=0;
    toadd->X[3]=0;
    #endif
    toadd->B[0]=0;
    toadd->B[1]=0;
    toadd->B[2]=0;
    toadd->B[3]=0;
    #ifdef THREE_DIMENSIONS
    toadd->B[4]=0;
    toadd->B[5]=0;
    #endif
    #ifdef FOUR_DIMENSIONS
    toadd->B[4]=0;
    toadd->B[5]=0;
    toadd->B[6]=0;
    toadd->B[7]=0;
    #endif
    //Add node at free position
    if(curr->left==NULL) {
        curr->left=toadd;
    }
    else if(curr->right==NULL) {
        curr->right=toadd;
    }
    else cout<<"Error, could not add node."<<j<<endl;
}

void set_trivial_param(saw_node* root)
{
int n;
saw_node* curr=root;
n=curr->n;
if(n>2){
curr->left->n=n/2;
curr->left->X[0]=n/2;
curr->right->n=n/2;
curr->right->X[0]=n/2;
set_trivial_param(curr->left);
set_trivial_param(curr->right);
}
return;
}

/*Set all node parameters to the initial configuration (straight rod configuration)*/
void set_initial_param(saw_node* root)
{
//Depth D, Level L, # of steps n.
saw_node* curr=root;
int D,L,n;
n=root->n;
D=ceil(log((double)n)/log(2.0));
for(L=1;L<D-1;L++) {
    if(n-pow(2,D-L)>=rint(ldexp(1,D-L-1))) {
    curr->left->X[0]=rint(ldexp(1,D-L));
    curr->right->X[0]=n-rint(ldexp(1,D-L));
    set_trivial_param(curr->left);
    curr=curr->right;
    }
    else {
    curr->left->X[0]=n-rint(ldexp(1,D-L-1));
    curr->right->X[0]=rint(ldexp(1,D-L-1));
    set_trivial_param(curr->right);
    curr=curr->left;
    }
}
}
/*Generates saw tree with n nodes.*/
struct saw_node *generate_saw_tree(int n)
{
int j;
if(n<=0){
cout<<"invalid saw-tree size:"<<n<<endl;
return NULL;
}
#ifdef TWO_DIMENSIONS
saw_node *root=new saw_node;
root->n=n;
root->X[0]=n; root->X[1]=0;
root->q=0;
//root->B[0]=0; root->B[1]=0; root->B[2]=n; root->[3]=0;
root->left=NULL;
root->right=NULL;
root->parent=NULL;
//Generate the n-1 internal nodes of the saw tree
for(j=1; j<n-1; j++) {
add_node(j,root);
}
set_initial_param(root);
return root;
#endif
}

/*Clear tree structure and free memory*/
void clear_saw_tree(struct saw_node *node)
{
	if(node != NULL) {
		clear_saw_tree(node->left);
		clear_saw_tree(node->right);
		delete node;
	}
}

/*Merge tree wl and wr to tree w*/
void merge(saw_node* wl, saw_node* wr, saw_node* w)
{
#ifdef TWO_DIMENSIONS
w->n=wl->n+wr->n;
#endif
}

/*-------------------*/
/*Internal Operations*/
/*-------------------*/

/*Symmetry operation applied to vector X.*/
void symX(int q,int* X,int* result)
{
#ifdef TWO_DIMENSIONS
//in total 8 symmetry operations
int temp;
switch(q) {
    case 0: result[0]=X[0]; result[1]=X[1]; break; //identity -> do nothing
    case 1: result[0]=X[1]; result[1]=X[0];  break; //reflect in 45deg
    case 2: result[0]=-X[1]; result[1]=-X[0]; break; //reflect in -45deg
    case 3: result[0]=X[1]; result[1]=-X[0]; break; //rotate by 90deg clockwise
    case 4: result[0]=-X[0]; result[1]=-X[1]; break; //rotate by 180deg clockwise
    case 5: result[0]=-X[1]; result[1]=X[0]; break; //rotate by 270deg clockwise
    case 6: result[1]=-X[1]; result[0]=X[0]; break; //reflect at x-axis
    case 7: result[0]=-X[0]; result[1]=X[1]; break; //reflect at y-axis
    default:
    cout<<"Wrong symmetry operation: "<<q<<endl;
}
#endif
}

void addX(int* point1, int* point2, int* result)
{
#ifdef TWO_DIMENSIONS
result[0]=point1[0]+point2[0];
result[1]=point1[1]+point2[1];
#endif
#ifdef THREE_DIMENSIONS
result[0]=point1[0]+point2[0];
result[1]=point1[1]+point2[1];
result[2]=point1[2]+point2[2];
#endif
#ifdef FOUR_DIMENSIONS
result[0]=point1[0]+point2[0];
result[1]=point1[1]+point2[1];
result[2]=point1[2]+point2[2];
result[3]=point1[3]+point2[3];
#endif
}

/*Symmetry operation q applied to bounding box B.*/
void symB(int q, int* B, int* result)
{
#ifdef TWO_DIMENSIONS
int temp;
switch(q) {
    case 0: result[0]=B[0]; result[1]=B[1]; result[2]=B[2]; result[3]=B[3]; break; //identity -> do nothing
    case 1: result[0]=B[1]; result[1]=B[0]; result[2]=B[3]; result[3]=B[2]; break; //reflect in 45deg
    case 2: result[0]=fmin(-B[1],-B[3]); result[1]=fmin(-B[0],-B[2]); result[2]=fmax(-B[1],-B[3]); result[3]=fmax(-B[0],-B[2]); break;//reflect in -45deg
    case 3: result[0]=B[1]; result[2]=B[3]; result[1]=fmin(-B[0],-B[2]); result[3]=fmax(-B[0],-B[2]); break; //rotate by 90deg clockwise
    case 4: result[0]=fmin(-B[0],-B[2]); result[2]=fmax(-B[0],-B[2]); result[1]=fmin(-B[1],-B[3]); result[3]=fmax(-B[1],-B[3]); break; //rotate by 180deg clockwise
    case 5: result[0]=fmin(-B[1],-B[3]); result[2]=fmax(-B[1],-B[3]); result[1]=B[0]; result[3]=B[2]; break; //rotate by 270deg clockwise
    case 6: result[0]=B[0]; result[2]=B[2]; result[1]=fmin(-B[1],-B[3]); result[3]=fmax(-B[1],-B[3]); break; //reflect at x-axis
    case 7: result[1]=B[1]; result[3]=B[3]; result[0]=fmin(-B[0],-B[2]); result[2]=fmax(-B[0],-B[2]); break; //reflect at y-axis
    default:
    cout<<"Wrong symmetry operation";
}
#endif
}

/*Translation by vector X applied to bounding box B.*/
void transB(int* X, int* B, int* result)
{
#ifdef TWO_DIMENSIONS
result[0]=B[0]+X[0]; result[1]=B[1]+X[0]; //x1-direction translation
result[2]=B[2]+X[1]; result[3]=B[3]+X[1]; //x2-direction translation
#endif
#ifdef THREE_DIMENSIONS
result[0]=B[0]+X[0]; result[1]=B[1]+X[0]; //x1-direction translation
result[2]=B[2]+X[1]; result[3]=B[3]+X[1]; //x2-direction translation
result[4]=B[4]+X[2]; result[5]=B[5]+X[2]; //x3-direction translation
#endif
#ifdef FOUR_DIMENSIONS
result[0]=B[0]+X[0]; result[1]=B[1]+X[0]; //x1-direction translation
result[2]=B[2]+X[1]; result[3]=B[3]+X[1]; //x2-direction translation
result[4]=B[4]+X[2]; result[5]=B[5]+X[2]; //x3-direction translation
result[6]=B[6]+X[3]; result[7]=B[7]+X[3]; //x4-direction translation
#endif
}

/* Fuses the boxes B1 and B2.*/
void fuseB(int* B1, int* B2, int* result)
{
#ifdef TWO_DIMENSIONS
result[0]=fmin(B1[0],B2[0]); result[1]=fmin(B1[1],B2[1]); // res_inf=inf(inf1,inf2)
result[2]=fmax(B1[2],B2[2]); result[3]=fmax(B1[3],B2[3]); // res_sup=sup(sup1,sup2)
#endif
#ifdef THREE_DIMENSIONS
result[0]=fmin(B1[0],B2[0]); result[1]=fmin(B1[1],B2[1]); result[2]=fmin(B1[2],B2[2]);// res_inf=inf(inf1,inf2)
result[3]=fmax(B1[3],B2[3]); result[4]=fmax(B1[4],B2[4]); result[5]=fmax(B1[5],B2[5]);// res_sup=sup(sup1,sup2)
#endif

#ifdef FOUR_DIMENSIONS
result[0]=fmin(B1[0],B2[0]); result[1]=fmin(B1[1],B2[1]); result[2]=fmin(B1[2],B2[2]); result[3]=fmin(B1[3],B2[3]);// res_inf=inf(inf1,inf2)
result[4]=fmax(B1[4],B2[4]); result[5]=fmax(B1[5],B2[5]); result[6]=fmax(B1[6],B2[6]); result[7]=fmax(B1[7],B2[7]);// res_sup=sup(sup1,sup2)
#endif
}


int main()
{
saw_node* root;
root=generate_saw_tree(17);
cout<<root->right->X[0]<<endl;
clear_saw_tree(root);
}
